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Abstract 

Efficient analysis and simulation of multiscale stochastic systems of chemical ki¬ 
netics is an ongoing area for research, and is the source of many theoretical and 
computational challenges. In this paper, we present a significant improvement 
to the constrained approach, which is a method for computing effective dynam¬ 
ics of slowly changing quantities in these systems, but which does not rely on 
the quasi-steady-state assumption (QSSA). The QSSA can cause errors in the 
estimation of effective dynamics for systems where the difference in timescales 
between the “fast” and “slow” variables is not so pronounced. 

This new application of the constrained approach allows us to compute the 
effective generator of the slow variables, without the need for expensive stochas¬ 
tic simulations. This is achieved by finding the null space of the generator of 
the constrained system. For complex systems where this is not possible, or 
where the constrained subsystem is itself multiscale, the constrained approach 
can then be applied iteratively. This results in breaking the problem down into 
finding the solutions to many small eigenvalue problems, which can be efficiently 
solved using standard methods. 

Since this methodology does not rely on the quasi steady-state assumption, 
the effective dynamics that are approximated are highly accurate, and in the case 
of systems with only monomolecular reactions, are exact. We will demonstrate 
this with some numerics, and also use the effective generators to sample paths 
of the slow variables which are conditioned on their endpoints, a task which 
would be computationally intractable for the generator of the full system. 
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1. Introduction 

Understanding of the biochemical reactions that govern cell function and 
regulation is key to a whole range of biomedical and biological applications and 
understanding mathematical modelling of gene regulatory networks has been an 
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area of huge expansion over the last half century. Due to the low copy numbers 
of some chemical species within the cell, the random and sporadic nature of 
individual reactions can play a key part in the dynamics of the system, which 
cannot be well approximated by ODEs[T3]. Methods for the simulation of such 
a system, such as Gillespie’s stochastic simulation algorithm (SSA)JH]i (or the 
similar Bortz-Kalos-Lebowitz algorithm]^ specifically for Ising spin systems), 
have been around for some decades. Versions which are more computationally 
efficient have also been developed in the intermediate vears|171 [7]. 

Unfortunately, their application to many systems can be very computation¬ 
ally expensive, since the algorithms simulate every single reaction individually. 
If the system is multiscale, i.e. there are some reactions (fast reactions) which 
are happening many times on a timescale for which others (slow reactions) are 
unlikely to happen at all, then in order for us to understand the occurrences of 
the slow reactions, an unfeasible number of fast reactions must be simulated. 
This is the motivation for numerical methods which allow us to approximate 
the dynamics of the slowly changing quantities in the system, without the need 
for simulating all of the fast reactions. 

For systems which are assumed to be well-mixed, there are many different 
approaches and methods which have been developed. For example the r-leap 
method [20] speeds up the simulation by timestepping by an increment within 
which several reactions may occur. This can lead to problems when the copy 
numbers of one or more of the species approaches zero, and a number of different 
methods for overcoming this have been presented pTIl 12] . 

Several other methods are based on the quasi steady-state assumption (QSSA). 
This is the assumption that the fast variables converge in distribution in a time 
which is negligible in comparison with the rate of change of the slow variable. 
Through this assumption, a simple analysis of the fast subsystem yields an ap¬ 
proximation of the dynamics of the slow variables. This fast subsystem can 
be analysed in several ways, either through analysis and approximation[6], or 
through direct simulation of the fast subsvstem[ll]. 

Another approach is to approximate the system by a continuous state-space 
stochastic differential equation (SDE), through the chemical Langevin equation 
(CLE)[T2]. This system can then be simulated using numerical methods for 
SDEs. An alternative approach is to approximate only the slow variables by an 
SDE. The SDE parameters can be found using bursts of stochastic simulation 
of the system, initialised at a particular point on the slow state space [T5]. the 
so-called “equation-free” approach. This was further developed into the con¬ 
strained multiscale algorithm (CMA)|5|, which used a version of the SSA which 
also constrained the slow variables to a particular value. Using a similar ap¬ 
proach to [6|, the CMA can similarly be adapted so that approximations of the 
invariant distribution of this constrained system can be made without the need 
for expensive stochastic simulations HU]. However, depending on the system, as 
with the slow-scale SSA, these approximations may incur errors. Work on how 
to efficiently approximate the results of multiscale kinetic Monte Carlo problems 
is also being undertaken in many different applications such as Ising models and 
lattice gas models [24]. 
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Analysis of mathematical models of gene regulatory networks (GRNs) is 
important for a number of reasons. It can give us further insight into how im¬ 
portant biological processes within the cell, such as the circadian clock[33] or 
the cell cycle[531 work. In order for these models to be constructed, we need 
to observe how these systems work in the first place. Many of the observation 
techniques, such as the DNA microarray [37], are notoriously subject to a large 
amount of noise. Moreover, since the systems themselves are stochastic, the 
problem of identifying the structure of the network from this data is very diffi¬ 
cult. As such, the inverse problem of characterising a GRN from observations 
is a big challenge facing our communitv|21|. 

One popular approach to dealing with inverse problems, is to use a Bayesian 
framework. The Bayesian approach allows us to combine prior knowledge about 
the system, complex models and the observations in a mathematically rigorous 
wav|29j. In the context of GRNs, we only have noisy observations of the concen¬ 
trations of species at a set of discrete times. As such, we have a lot of missing 
information. This missing data can be added to the state space of quantities that 
we wish to infer from the data that we do have. This complex probability distri¬ 
bution on both the true trajectories of the chemical concentrations, and on the 
network itself, can be sampled from using Markov chain Monte Garlo (MCMC) 
methods, in particular a Gibb’s sampler [IB]. Within this Gibb’s sampler, we 
need a method for sampling a continuous path for the chemical concentrations 
given a guess at the reaction parameters, and our noisy measurements. Exact 
methods for sampling paths conditioned on their endpoints have been developed 

HSlEll. 

In other applications, methods for path analysis and path sampling have 
been developed, for example discrete path sampling databases for discrete time 
Markov chains|33], or where the probability of paths, rather than that of trajec¬ 
tories of discrete Markov processes can be used to analyse behaviour |30]. In [T3| . 
a method for transition path sampling is presented for protein folding, where 
the Markov chain has absorbing states. Other approaches for coarse-graining 
transition path sampling in protein folding also exist [3|. Other methods also ex¬ 
ist for the simulation of rare events where we wish to sample paths transitioning 
from one stable region to another]!]. 

The problems become even more difficult when, as is often the case, the 
systems in question are also multiscale. This means that these inverse problems 
require a degree of knowledge from a large number of areas of mathematics. 
Even though many of the approaches that are being developed are currently 
out of reach in terms of our current computational capacity, this capacity is 
continually improving. In this paper we aim to progress this methodology in a 
couple of areas. 

1.1. Conditioned path sampling methods 

We will briefly review the method presented in m for the exact sampling 
of conditioned paths in stochastic chemical networks. Suppose that we have a 
Markov jump process, possibly constructed from such a network, with a gener¬ 
ator Q. The generator of such a process is the operator Q such that the master 
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[1] Define a dominating process to have transition rates given by the matrix 

+ 

[2] This process has uniformly distributed reaction events on the time interval 

The number r of such events is given by Q. 

[3] Once r = f has been sampled, the type of each event must be decided, 
by sampling from the distribution ([^, starting with the first event. An 
event which corresponds to rate rrii^i indicates that no reaction event has 
occurred at this event. 

[4] Once all event types have been sampled, we have formed a sample from 
the conditioned path space. 


Table 1: A summary of the methodology presented in m, for sampling paths of Markov- 
modulated Poisson processes, conditioned on their endpoints. 


equation of the system can be expressed as 


dp 

dt 


= Gp, 


where p is the (often infinite dimensional) vector of probabilities of being in 
a particular state in the system. We wish to sample a path, conditioned on 
X(to) = xo and X(ti) = xi. Such a path can be found by creating a domi¬ 
nating process (i.e. a process whose rate is greater than the fastest rate of any 
transitions of the original system) with a uniform rate. 

We define the rate to be greater than the fastest rate of the process with 
generator Q, so that 

p > max0i_i. 
i 

Then we define the transition operator of the dominant process by: 

M = -g + i. 
p 


We can then derive the number of reaction events Nu of the dominating process 
in the time interval by: 


n^u = r) 


eyip{-pt){ptY 

[exp(^t)],ro,xt 


( 1 ) 


Here the notation [-Ja.b denotes the entry in the matrix with coordinates (a, b) G 
N^. A sample is taken from this distribution. The times {t*, ^ 2 ) ■ • ■ t*} of all of 
the r reaction events can then be sampled uniformly from the interval [tQ,ti\. 
The only thing that then remains is to ascertain which reaction has occurred at 
each reaction event. This can be found by computing, starting with A (to) = xq, 
the probability distribution defined by: 


¥{X{t*)) = x\X{t*_,) = x*_„X{tY = xY 


Yx,Xi 


( 2 ) 
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This method, summarised in Table exactly samples from the desired distri¬ 
bution, but depending on the size and sparsity of the operator C/, it can also 
be very expensive. In the context of multiscale systems with a large number of 
possible states of the variables, the method quickly becomes computationally 
intractable. 

1.2. Summary of Paper 

In Section we introduce a version of the Constrained Multiscale Algorithm 
(CMA), which allows us to approximate the effective generator of the slow pro¬ 
cesses within a multiscale system. In particular, we explore how stochastic 
simulations are not required in order to compute a highly accurate effective 
generator. In Sectionwe consider the differences between the constrained ap¬ 
proach, and the more commonly used quasi-steady state assumption (QSSA). 
In Section we describe how the constrained approach can be extended in an 
iterative nested structure for systems for whose constrained subsystem is itself 
a large intractable multiscale system. By applying the methodology in turn to 
the constrained systems arising from the constrained approach, we can make 
the analysis of highly complex and high dimensional systems computationally 
tractable. In Sectionj^ we present some analytical and numerical results, aimed 
at presenting the advantages of the CMA over other approaches. This includes 
some examples of conditioned path sampling using effective generators approx¬ 
imated using the CMA. Finally, we will summarise our findings in Section 

2. The Constrained Multiscale Algorithm 

The Constrained Multiscale Algorithm was originally designed as a mul¬ 
tiscale method which allowed us to compute the effective drift and diffusion 
parameters of a diffusion approximation of the slow variables in a multiscale 
stochastic chemical network. The idea was simply to constrain the original dy¬ 
namics to a particular value of the slow variable. This can be done through a 
simple alteration of the original SSA by Gillespie [TH]. First, a (not necessarily 
orthogonal) basis is found for the system in terms of “slow” and “fast” vari¬ 
ables, [S = [S'!, 5'2,...], F = [Fi, F 2 ,...]]. Slow variables are not affected by the 
most frequently firing reactions in the system. Then, as shown in [9], the SSA 
is computed as normal, until one of the slow reactions (a reaction which alters 
the value of the slow variable(s)) occurs. After the reaction has occurred, the 
slow variable is then reset to its original value, in such a way that the fast vari¬ 
ables are not affected. This is equivalent to projecting the state of the system, 
after each reaction, back to the desired value of the slow variable, whilst also 
preserving the value(s) of the fast variable(s). The constrained SSA is given in 
Table Here the ai(X(t)) denote the propensity of the reaction Ri when the 
system is in state X(t) = [Xi(t), A 2 (t),...], where Atoi(X(t) is the probability 
that this reaction will fire in the infinitesimally small time interval (t, t + At) 
with 1 ^ At > 0. The stoichiometric vectors i^i denote the change in the state 
vector X(t) due to reaction Ri firing. 
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In order to describe the constrained approach, we first introduce some defi¬ 
nitions that will be helpful. 

Definition 2.1. Constrained Projector: Given a basis of the state space X = 
[Xi,X 2 ,... ,Xn] with Nf fast variables F = [Fi, Fj,...,F/v, 
ables S = S' 2 ,. •., Sn^], the constrained projeetor Vs ■ Nq —>■ for a given 

value of S preserves the values of the fast variables, whilst mapping the values 
of the slow variables to S .■ 

Fs([S,F]) = [S,F] V([S,F])e<. (3) 

Definition 2.2. Constrained Stoichiometric Projector: Given a basis of the 
state space X = [Xi, X 2 ,..., Xn] with Nf fast variables F = [Fi, F 2 ,..., FjVj] 
and Ns slow variables S = [Fi, S' 2 ,..., the constrained stoichiometric pro¬ 

jector V : Nq —>■ Nq maps any non-zero elements of the slow coordinates to 
zero, whilst preserving the values of the fast coordinates: 

F([S,F]) = [0,F] V([S,F])gN^. (4) 


and Ns slow vari- 


Definition 2.3. Constrained Subsystem: Given a system with Nji reactions 
Ri, R 2 ,. ■ ■, Rnr with propensity functions ai(X) and stoichiometric vectors 
Vi G N^, the constrained subsystem is the system that arises from applying 
the constrained projector Vs to the state vector after each reaction in the sys¬ 
tem. This is equivalent to applying the constrained stoichiometric projector V 
to each of the stoichiometric vectors in the system. This may leave some reac¬ 
tions with a null stoichiometric vector, and so these reactions can be removed 
from the system. This projection can lead to aphysical systems where one or 
more variables may become negative; in these cases we set the propensities of 
the offending reactions at states where a move to a negative rate is possible, to 
zero. 

Let us illustrate this using an example which we shall be using later in the 
paper. 


Ri : 

0 - 

^ Ai, 


R 2 : 

X 2 

^ 0 , 

(5) 

F 3 : 

Xi 

-^2, 


! 

X 2 

-S Ai. 



In certain parameter regimes, this system is multiscale, with reactions F 3 and 
F 4 occurring many times on a time scale for which reactions i?i and R 2 are 
unlikely to happen at all. The variable S' = Xi -|- X 2 is unaffected by these 
fast reactions, and as such is a good candidate for the slow variable which we 
wish to analyse. A discussion about how the fast and slow variables could be 
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[1] Define a basis of the state space in terms of slow and fast variables. 

[2] Initialise the value of the state, X(to) = x. 

[3] Calculate propensity functions at the current state Q!i(X(t)). 

[4] Sample the waiting time to the next reaction in the system 


T 


log (u) 
ao(X(t))’ 


M 

where ao{X{t)) = ^ ak{X{t)), 

k=l 


u~t/([0,l]). 


[5] Choose one j G {1,..., M}, with probability aj/ao, and perform reaction 
Rj, with stoichiometry which has been projected using the constrained 
stoichiometric projector: 

x{t + T) = x{t) + r{i^,). 

[ 6 ] Repeat from step [3], 


Table 2: The Constrained Stochastic Simulation Algorithm (CSSA) using the constrained 
stoichiometric projector given in Definition \2.2\ Simulation starts with S = s where s is a 
given value of the slow variable. 


identified is given in Section]^ We have two choices for the fast variable, either 
F = Xi or F = X 2 , in order to form a basis of the state space along with the 
slow variable S. As detailed in [S], it is preferable (although not essential) to 
pick fast variables that are not involved in zeroth order reactions. Therefore, 
in this case, we choose F = X 2 . Following the projection of the stoichiometric 
vectors using the constrained projector, the constrained system can be written 
in the following way: 


Cl : 

Xi 

+ X 2 — s, 


i?2 : 

X 2 

-^ -VI, 


i?3 : 

Xi 

-5" -V2, 

(6) 

i?4 : 

X 2 

A Ai. 



Note that reaction Ri has disappeared completely, since only involves changes 
to the slow variable, and as such after projection, the stoichiometric vector is 
null, and the reaction can be removed. The stoichiometry of reaction R 2 has 
been altered as it involves a change in the slow variables. If this reaction occurs, 
the slow variable is reduced by one. We are not permitted to change the fast 
variable X 2 in order to reset the slow variable to its original value, and therefore 
we must increase Xi by one, giving us a new stoichiometry for this reaction. 

In the original CMA, statistics were taken regarding the frequency of the 
slow reactions, at each point of the slow domain, and were used to construct 
the effective drift and diffusion parameters of an effective diffusion |9l|8] process. 
However, this constrained approach can also be used to compute an effective 
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generator for the discrete slow process, as we will now demonstrate. The CMA 
can be very costly, due to the large computational burden of the stochastic 
simulations of the constrained system. In this section, we will introduce a 
method for avoiding the need for these simulations, whilst also significantly 
improving accuracy. 

The constrained systems can often have a very small state space (which 
we will denote r(s)), since they are constrained to a single value of the slow 
variables. For example, for the constrained system Q, there are only 
possible states. Such a system can easily be fully analysed. For example, the 
invariant distribution can be found by characterising the one-dimensional null 
space of the generator matrix of the constrained process. For small to medium¬ 
sized systems, this is far more efficient than exhaustive Monte Carlo simulations. 

For other systems with larger constrained state spaces, stochastic simulation 
may still be the best option, although in Section]^ we show how the constrained 
approach can be applied iteratively until the constrained subsystem is easily 
analysed. 

Suppose that we have a constrained system with Nf fast variables, Fi, F 2 ,..., Fn^ . 
The generator for the constrained system with S = s is given by Gpis). Since the 
system is ergodic, there is a one-dimensional null space for this generator. This 
can be found by using standard methods for identifying eigenvectors, by search¬ 
ing for the eigenvector corresponding to the eigenvalue equal to zero. Krylov 
subspace methods allow us to find these eigenvectors with very few iterations. 
Suppose we have found such a vector v = [ui, W 2 ) ■ • •]> such that 

Gf{s)v = 0. 

Then our approximation to the invariant distribution of this system is given by 
the discrete probability distribution represented by the vector 

p(s) = [pi(s),P 2 (s),...] = 

Our aim is now to use this distribution to find the effective propensities of the 
slow reactions of the original system. 

Suppose that we have Ng slow reactions in the original system. Each has an 
associated propensity function ai(S', F), 02 ( 5 , F),..., (5', F). We now sim¬ 

ply want to find the expectation of each of these propensity functions with 
respect to the probability distribution p(s): 

E{a^{S,-)) = Y^Pi{s)ai{SJ). (7) 

i 

Having computed this expectation for all of the slow propensities, over all re¬ 
quired values of the slow variable, then an effective generator for the slow vari¬ 
able can be constructed. 

3. Comparing the CMA and QSSA approaches 

A very common approach to approximating the dynamics of slowly changing 
quantities in multiscale systems, is to invoke the quasi steady-state assumption 



[1] For each value of the slow variable S' = s S 17, compute the generator Qs 
of the constrained subsystem. 

[2] Find the zero eigenvector v = [vi,V 2 , ...] of t/s, and let p(s) = 

[3] Approximate the effective propensities at each point s S 17 using Q. 

[4] Construct an effective generator Q of the slow processes of the system 
using these effective propensities. 


Table 3: The CMA approach to approximating the effective generator Q of the slow variables 
on the (possibly truncated) domain S £0., without the need for stochastic simulations. 


(QSSA). The assumption is that the fast and slow variables are operating on 
sufficiently different time scales that it can be assumed that the fast subsystem 
enters equilibrium instantaneously following a change in the slow variables, and 
therefore is unaffected by the slow reactions. This assumption means that if the 
fast subsystem’s invariant distribution can be found (or approximated), then 
the effective propensities of the slow reactions can be computed. However, as 
demonstrated in [S], this assumption incurs an error, and for systems which do 
not have a large difference in time scales between the fast and slow variables, 
this error can be significant. 

The CMA does not rely on the QSSA, and is able to take into account 
the effect that the slow reactions have on the invariant distribution of the fast 
variables, conditioned on a value of the slow variables. In a true fast-slow 
system, this will yield the same results as the QSSA, but for most systems of 
interest, the constrained approach will have a significant increase in accuracy. 
If we follow the approach outlined in Table we don’t even need to conduct 
any stochastic simulations to approximate the effective dynamics. 

The assumptions for the CMA are weaker than the QSSA, namely that 
we assume that the dynamics of the slow variable(s) can be approximated by 
a Markov-modulated Poisson process, independently of the value of the fast 
variables. This means that we have made the assumption that the current value 
of the fast variables has no effect on the transition rates of the slow variables 
once a slow reaction has occurred. This is subtly weaker than the QSSA, and 
importantly the effect of the slow reactions on the invariant distribution of the 
fast variables is accounted for. Note that this may necessitate a slow variable 
which has more than one dimension, for example in oscillating systems for which 
the effective dynamics cannot be approximated by a one dimensional Markov 
process. Consideration of such systems is an area for future work. 

4. The Nested CMA 

There will be many systems for which the constrained subsystem is itself a 
highly complex and multiscale system. In this event, it will not be feasible to find 
the null space of a sensibly truncated generator for the constrained subsystem. 
Therefore, we need to consider how we might go about approximating this. 
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Fortunately, we already have the tools to do this, since we can iteratively apply 
the CMA methodology to this subsystem. This is analogous to the nested 
strategy proposed in the QSSA-based nested SSA[TT]. 

This nested approach allows us to reduce much more complex systems in 
an accurate, computationally tractable way. The problem of finding the null 
space of the first constrained subsystem is divided into hnding the null space of 
many small generators, through further constraining. An example of this nested 


approach will be demonstrated in Section 5.3 


5. Examples 

In this section we will present some analytical and numerical results produced 
using the CMA approach for three different examples. In order to give an 
indication of the computational cost of the algorithms, we include the runtime of 
certain operations. All numerics were performed using MATLAB on a mid-2014 
MacBook Pro. Disclaimer: the implementations used are not highly optimised, 
and these runtimes are purely given as an indication of the true costs of a well 
implemented version. 


5.1. A Simple Linear System 

First we consider a simple linear system, in order to demonstrate that the 
CMA approximation of the effective generator of the slow variable is exact in 
the case of systems with only monomolecular reactions, which is in contrast to 
the approximation found using a more standard QSSA-based approach. Let us 
illustrate this by returning to the example given by the linear system ([^, first 
analysing it using the QSSA. 


Ri : 

0 - 

fel, V 
— Ai, 

i ?2 : 

^2 

^ 0, 

R3 : 


-^ A2 

R 4 : 

^2 

A Ail 


We will consider this system in the following parameter setting: 

kiV = 20, fca = 1, fca = 5, ki = 5. (8) 

Here V denotes the volume of the well-mixed thermally-equilibrated reactor. 


5.1.1. QSSA-based analysis 

The QSSA tells us that the fast subsystem (made up of reactions and R^) 
reaches probabilistic equilibrium on a timescale which is negligible in comparison 
with the timescale on which the slow reactions are occurring. Therefore we may 
treat this subsystem in isolation with fixed S: 


Xi 



S — Xi X2. 
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This is a very simple autocatalytic reaction system, for which a great deal 
of analytical results are available. For instance, we can compute the invariant 
distribution for this system [22]. which gives us that X 2 is a binomial random 
variable 


^2 


B 



Therefore, we can compute the conditional expectation E(X2|5') = in this 

fast subsystem, and use this to approximate the effective rate of reaction i? 2 - 
Therefore, the effective slow system is given by the reactions: 


0 



S 



(9) 


where 


ki = ki, 


A:2E(^2) ^2^3 

S k3 + ki 


We can compute the invariant distribution for this effective system [52], which 
in this instance is a Poisson distribution: 


S 


/ kiVik3 + ki) \ 
1 k2k3 ) 


( 10 ) 


We can quantify the error we have made in using the quasi-steady state as¬ 
sumption by, for example, comparing this distribution with the true invariant 
distribution. Once again, using the results of |22j , we can compute the invariant 
distribution of the system ([^, which is a multivariate Poisson distribution: 

[Xl,X2] ~ 7 ^(Ai,A2), 

where Ai = , and A 2 = Trivially one can compute the marginal 

distribution on the slow variable S: 


^ \n \s—n 

P(S=.) = 


■ exp(-(Ai -h A 2 )), 


n! (s — n)! 

exp(-(Ai -h A 2 )). 


n—0 

(Ai -I- A2) 


s! 

Therefore S is also a Poisson variable with intensity A = A 1 -I-A 2 = ’ 


which differs from the intensity approximated invariant density (101 by 
Note that k^ is one of the fast rates, and kiV is one of the slow rates, and 
therefore as the difference in timescales of the fast and slow reactions increases, 
this error decreases to zero, so that the QSSA gives us an asymptotically exact 
approximation of the slow dynamics. For systems with a finite timescale gap, 
the QSSA approximation will incur error over and above the error incurred in 
any approximation of the marginalised slow process by a Markov process. 
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5.1.2. CM A analysis 

For comparison, let us compute approximations of the effective slow rates 
by using the CMA. The CMA for this system tells us that we need to analyse 
the constrained system §)• 


Cl : 

^1 

+ X 2 

= 5 

i?2 : 

^2 


^1, 

i?3 : 

Xi 


^2, 

i?4 : 

X 2 

fc4 ^ 

Xi. 


The constrained system in this example only contains monomolecular reactions, 
and as such can be analysed using the results of [22] . The invariant distribution 
for this system is a binomial, such that 




^2 + ^3 + ^4 

Using this, we can compute the effective propensity of reaction R 2 , 

k2hS 


02(5) =fc2E(X2|5) = 


^^2 T ^3 “ t “ ^^4 


giving us the effective rate k 2 = The invariant distribution of 

with this effective rate for k 2 is once again a Poisson distribution with intensity 


A = 


kiV{k2 


ki) 




2«'3 


which is identical to the intensity of the true distribution on the slow vari¬ 
ables. In other words, for this example, the CMA produces an approximation 
of the effective dynamics of the slow variables for this system, whose invari¬ 
ant distribution is identical to the marginal invariant distribution of the slow 
variables in the full system. The constrained approach corrects for the effect 
of the slow reactions on the invariant distribution of the fast variables. In this 
and other examples of systems with monomolecular reactions, the constrained 
approach gives us a system whose invariant distribution is exactly equal to the 
marginal distribution on the slow variables for the full system. Another example 
is presented in Section for which the constrained system is itself too large to 
easily compute expectations directly through its generator, and requires another 
iteration of the CMA to be applied. 

For this example, we did not even need to compute the invariant distri¬ 
butions of the constrained systems numerically. In Section 5.2 we will come 


across a system for which it is necessary to numerically compute the invariant 
distribution of the constrained system. 
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Figure 1: Approximations of the invariant distribution of the slow variable S = Xi + X 2 of 
system with parameters (|^ through marginalisation of the distribution of the full system 
(histogram), of the effective generator computed using the CM A (solid line) and computed 
using the QSSA (dashed line). 



5.1.3. Comparison of approximation of invariant densities 

Figure [ 1 ] shows the invariant distributions of the slow variables S = Xi +X 2 
in the parameter regime ([^, computed by marginalising the invariant distribu¬ 
tion of the full system, and from the CMA and QSSA as outlined above. The 
CMA exactly matches the true distribution, as both are Poisson distributions 
with rate A = 44. The QSSA incorrectly approximates the effective rate of i? 4 , 
and as such is a Poisson distribution with rate A = 40. The relative error of the 
CMA for this problem is zero, and for the QSSA is 4.322 x 10“^. 

5.1.4- Conditioned path sampling using effective generators 

The approaches described in Section o hit problems when the system for 
which we are trying to generate a conditioned path is multiscale. In a multiscale 
system, the rate p of the dominating process will be very large, and as such 
the number of reaction events will be large, even if the path we are trying to 
sample is short. Therefore AP" is likely to be a full matrix, and the number of 
calculations of ([^ will be large. Moreover, the size of M is also likely to be 
large, since for each value S' = s of the slow variable, there are many states, 
one for each possible value of the fast variable. All of these factors make the 
problem of computing a conditioned path in such a scenario computationally 
intractable. 

Considering once more the system ([^, naturally we cannot store the actual 
generator of this system, since the system is open and as such the generator 
is an infinite dimensional operator. However, the state space can be truncated 
carefully in such a way that the vast majority of the states with non-negligible in¬ 
variant density are included, but an infinite number of highly unlikely states are 
presumed to have probability zero. Note that this means that we are effectively 
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sampling paths satisfying S{to) = si, S'(ti) = S 2 conditioned on S(t) G Vt. 
However, even with careful truncation the number of states can be prohibitively 
large. 

Suppose that we truncate the domain for this system to 
n = {[Xi,X2]\XuX2 G {0,1,...,200}. 

This truncated system has 201^ = 40401 different states, and therefore the gen¬ 
erator Q G k40401x40401^ Although this matrix is sparse, the matrix exponential 
required in Q is full, as is M’' for moderate r € N. A full matrix of this size 
stored at double precision would require over 13GB of memory. So even for this 
system, the most simple multiscale system that one could consider, the problem 
of sampling conditioned paths is computationally intractable. 

In comparison, suppose that we use a multiscale method such as the CMA to 
approximate the effective rates of the slow reactions. Then, for the same H, we 
only have 401 possible states of the slow variable, a reduction of 99.25%. The 
effective generator Q G k401x401 then only require 1.29MB to be stored 

as a full matrix in double precision. The dominating process for this system 
must now have rate p > 201.4, instead of p > 1220, which is over 6 times bigger. 
This means far fewer calculations of (§. What is more, as we saw in Section 
|5.1.2[ for some systems the CMA exactly computes the effective dynamics of 
the slow variables, with no errors. 

The system ([^, in order to highlight more effectively the differences between 
the CMA and a QSSA-based approach, is in a parameter range Q, for which 
the difference in time scales between the “fast” and “slow” variables is relatively 
small, and of course for systems with larger timescale difference, the difference 
in p between the full and effective generators would be far larger. 

Naturally, this approach only allows us to sample the paths of the slow 
variables. However, the fast variables, if required, can easily be sampled after 
the fact, using an adapted Gillespie approach which samples the fast variables 
given a trajectory of the slow variables. 

As we have just demonstrated in the previous section, the CMA can be used 
to compute an effective generator for the slow variable S = X 1 +X 2 in the system 
([^ , with parameters (|^ , whose invariant distribution is exactly that of the slow 
variable in the full system without the multiscale reduction. Moreover, this can 
be achieved with no Monte Carlo simulations, since the constrained subsystem 
contains only monomolecular reactions, and as such its invariant distribution 
can be exactly computed [52]. 

At this juncture, we simply need to apply the method of Fearnhead and 
Sherlock[T6] to the computed effective generator in order to be able sample paths 
conditioned on their endpoints. Suppose we wish to sample paths conditioned 
on S{to = 0) = 44 = = 10). The invariant distribution of this system, 

as shown previously in this paper, is a Poisson distribution with mean A = 
~ Therefore, we are attempting to sample paths which start 
and finish at the the mean of the invariant distribution, which in itself is not a 
particularly interesting thing to do, but it will allow us to highlight again the 
advantages of using the CMA over QSSA-based approaches. 
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Since the system is open, we are required to truncate the domain in order 
to be able to store and manipulate the effective generator. We truncate the 
domain to = {[Xi,X2]|5' = Xi + X 2 < 400}. Therefore we aim to sample 
paths 

{S{t), t G [0,10] I 5(0) = 44 = 5(10), S{t) G flVt G [0,10]}. 

As the number of possible states of the slow variable is relatively small, it 
was possible to compute and store full matrices for M’' as required in Q and 
([^ for r G 1,2,..., 2369. r has an upper bound of 2369 as the cumulative mass 
function for the probability distribution Q is within machine precision of one 
at r = 2369. Storing all powers of the matrices is clearly not the most efficient 
way to implement this algorithm, but for this example was possible without any 
intensive computations, and with minimal numerical error. We will present a 
more efficient approach in the next section. 




(a) (b) 

Figure 2: (a) 10 trajectories of the slow variable conditioned on S'(O) = 44 = 5(10), sam¬ 
pled using the CM A approximate effective generator, (b) Mean and standard deviation of 
1000 trajectories of the slow variable conditioned on 5(0) = 44 = 5(10), sampled using the 
approximate effective generator from both the QSSA (blue plots) and CMA (red plots). 


Figure 5.1.4| (a) shows 10 example trajectories sampled using the the condi¬ 
tioned path sampling algorithm with the CMA approximation of the effective 
generator of the slow variable. We also implemented exactly the same approach 
using the QSSA approximation of the effective generator. The mean and stan¬ 
dard deviation of 1000 paths sampled using both methods is plotted in Figure 
5.1.4 (b). Since the paths are conditioned to start and finish at the mean of 


the system’s monomodal invariant distribution, we would expect the mean to 
converge to a constant 5 = 44 as we sample more paths. 

This appears to be the case for the paths sampled using the CMA effective 
generator, which is what we would hope since this generator preserves the true 
mean of the slow variables, as demonstrated in the previous section. 

The QSSA, as has also been demonstrated in Section [5.1.1[ does not correctly 
preserve the invariant distribution of the slow variables, and underestimates the 
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mean value of the invariant distribution. This can be seen in 5.1.4 (right), where 
the mean value of the path dips in the middle of the trajectory as it reverts to the 
mean of the invariant distribution of the QSSA approximation, before increasing 
towards the end of the trajectory in order to satisfy the condition S = 44. 

This demonstrates that the accuracy of the approximation of the effective 
dynamics has a knock-on impact, as one would expect, to the accuracy of the 
conditioned path sampling. It would be preferable, naturally, if we could com¬ 
pare path statistics of the multiscale approaches to that of conditioned paths 
statistics of the full system. However, this is simply not feasible, due to the 
size of the matrices, even for the truncated domain Q. Instead, this does suc¬ 
ceed in demonstrating that these methods make conditional path sampling of 
the slow variables a possibility, where it was computationally intractable pre¬ 
viously. Since the rates could be explicitly calculated for this simple example, 
the effective generators could be produced in the order of 10“ 3 seconds for the 
domain S' G {0,1,..., 400}. The set up process for the path sampling, involving 
finding the probabilities in 0 and computing the required powers of A4 took 
around 90 seconds. After this, each path took a third of a second to sample. 


5.2. A Bistable Example 

Sampling of paths conditioned on their endpoints is an integral part of some 
approaches to Bayesian inversion of biochemical data. A Gibb’s sampler can be 
used to alternately update the network structure and system parameters, and 
the missing data (i.e. the full trajectory), sampled for example using the method 
found in m- However, efficient methods to sample paths of multiscale systems 
may also be useful in other areas. For instance, it may allow us to sample paths 
which make rare excursions, or large deviations from mean behaviour. This 
forms part of the motivation for considering the next example. 

Let us consider the following chemical system, which in certain parameter 
regimes exhibits bistable behaviour. 


, i?2 : 

A 2 ^ . 

k2 


ks 

i?3, i?4 : 

0 ^ A; 

^4 

i?5, i?6 : 

Ai+Ai 

i?7 : 

A 2 0 


X1+X2 


( 11 ) 


A, 


fee 


In particular, we consider parameter regimes where the occurrence of reactions 
i ?5 and i?6 are on a relatively faster timescale than the other reactions. The 
following is just such a parameter regime: 


ki = 142, 
/c4 = 92.8, 


^ = 1 
H ’ 


kr. 


ksV = 880, 
fcg = 500, kr = 6. 


( 12 ) 
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As before, V denotes the volume of the well-mixed thermally-equilibrated reac¬ 
tor. 



X, 


(a) (b) 

Figure 3: faj A log plot of an approximation ttq of the invariant distribution on the slow 
variable S = X\+2X2 of system 03 with parameters ( |5.2[ |, demonstrating the bistable nature 
of the system. Approximation was computed by finding the null space of the full generator of 
the system on the truncated domain {0,1, •. ., 800} x {0,1,.. •, 1200}. (b) Proportion of total 
propensity (Xi, X 2 ) attributed to the fast reactions and Rq^ given by 


That said, this parameter regime is one in which the use of the QSSA will 
create significant errors, since the timescale gap is not very large in all parts of 
the domain as demonstrated in Figure Figure (a) shows a highly accurate 
approximation of the invariant distribution of the full system, found by com¬ 
puting the null space of the full generator for the system truncated to the finite 
domain O = {(xi, 0 : 2 ) G {0,1,..., 800} x {0,1,..., 1200}. The zero eigenvec¬ 
tor of this truncated generator was found using standard eigenproblem solvers, 
then normalised. Since this system has 2nd order reactions, its invariant den¬ 
sity cannot in general be written in closed form, and as such, we could use this 
approximation on the truncated domain in order to quantify the accuracy of the 
multiscale approaches. This plot demonstrates the bistable nature of this sys¬ 
tem, which can take a long time to switch between the two favourable regions. 
This example has been chosen in order that such an approximation can still be 
computed in order to check the accuracy of the approach. 

Figure]^ (b) shows the proportion of the total propensity for each state which 
is attributed to the fast reactions, and i? 5 , given by: 


Pn,,RAXi,X2) 


Q:5(^1! Ar2) -I- aejXi, X 2 ) 
ao{Xi,X2) 


o^sjXi, X2) + aejXi, X2) 
j:Z,adX„X2) 


(13) 

This proportion, which is a measure of the gap in timescales between the “fast” 
reactions R 5 and Re, and the rest of the reactions, varies across the domain. 
We can approximate the expected proportion of propensity attributed to the 
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fast reactions: 


^iPR5,Re) — ^ PR5,Rei^l^^2)'^QiXi,X2), 

(Xi,X2)en 


where ttq is the approximate invariant density of the full generator on the trun¬ 
cated domain il. In this system with parameters (5.2), E(P/{ 5 _/{g) = 0.6941, i.e. 
the expected proportion of all reactions which are either of type or i?g is 
69.41%. As such, although reactions i ?5 and Rq are occurring more frequently 
than other reactions, there is not a stark difference in timescales, as we might 
expect in a system for which the QSSA yields a good approximation. The “fast” 
reactions in this example are reactions i ?5 and Rq, and as such, S = Xi + 2 X 2 
is a good choice of slow variable, since this quantity is invariant to these fast 
reactions. 


5.2.1. The QSSA Approach 

By applying the QSSA to the system ( 0 , we can approximate the effective 
rates of the slow variables by considering the fast reactions in isolation. The 
fast subsystem is given by the reactions P 5 and Rq: 

Cl : Xi+ 2 X 2 = S, (14) 

fes 

i?5,i?6 : Ai -I-Ai ^^ X 2 . 

fee 

Lines denoted by Ci in this and what follows denotes a constraint. It is impor¬ 
tant to keep a track of these constraints, since each one reduces the dimension 
of the state space by one. 

For a fixed value of ^ = Ai -|- 2 A 2 S {0,1,..., 5'max}, we wish to hnd the 
generator for the process X 2 (or equivalently Ai = 5" — 2 A 2 ) within this fast 
subsystem. The generator can be found by considering the master equation for 
each state X 2 = i: 

^ = -{a^i^S - 2i,i) + ae{S - 2i,i))pi +a^^S - 2{i-l),i - l)pi_i 

at 

+ aeiS-2{i + l),i)pi+i, 


where Pi{t) is the probability of X2{t) = i. Putting this set of differential 
equations into vector form gives us: 


dP 

dt 


GP, 


where Q is the generator of the fast subsystem (14). Note that since we are 
restricted to states such that Xi + X 2 = S for some value of S', there are only 
J possible different states, and as such 1 / £ J ^ LvJ, Even for moderately 
large values of S, the one-dimensional null space of such a sparse matrix is not 
computationally expensive to find, and when normalised gives us the invariant 
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density of X 2 (and therefore Xi if required). This invariant density can then 
be used to compute the expectation of the propensities of the slow reactions of 
the system for the state S' as in Q, and in turn be entered into the (truncated) 
effective generator for the slow variable. 


5.2.2. The Constrained Approach 

When using the CMA, the methodology is largely the same as was described 
for the QSSA-based approach in the last section. The only real difference lies in 
the subsystem which is analysed in order to compute the invariant distribution 
of the fast variables conditioned on the value of the slow variable. As we have 
done previously, we will consider each of the reactions in the system in turn, 
constraining the value of the slow variable to a particular value, whilst being 
sure not to change the value of the fast variables. There are two choices for 
the fast variable, in order to form a basis of the state space along with the slow 
variable S, but as explained in detail in [9], F = A 2 is the best choice, since 
there is a zeroth order reaction involving Xi, which can lead to an unphysical 
constrained subsystem, if this is chosen as the fast variable. 

With this choice of fast variables, the hrst four reactions all disappear in 
the constrained subsystem. This is because none of these reactions alter the 
fast variable, and as such the constrained stoichiometric projector maps their 
stoichiometric vectors to zero, and therefore reactions Ri, R2, R3, R4 have no 
net effect on the constrained subsystem. 

Reaction i ?7 differs in that it causes a change in the fast variable X 2 . The 
projector in this case maps the stoichiometric vector to [— 2 , 1 ]^ and therefore 

the net effect of reaction R 7 is equivalent to X 2 —A Xi + Xi. This leads to 
the following constrained system: 

Cl : Ai + 2 A 2 = S, 

fcs 

R5,i?6 : Xi + Xi ^ A 2 , 

fee 

i?7 : X2^Xi+Xi. 


Note that since reactions Rq and Rj have the same stoichiometry, this system 
can be simplihed by removing R 7 and adding its rate to Rq: 


Cl : 

Xi + 2 A 2 — S, 


i?5, i?6 : 

ka 

Xi+Ai ^ A 2 . 
fce+fc? 

(15) 


For every fixed value oi S G {0,1,..., Fmax}, the generator for (15) can be 
found following the same approach as in the previous section, the only difference 
being the altered rate for reaction Rg • Following this methodology, an effective 
generator Q can be computed. 


19 




s = + 2 X 2 


Figure 4: Approximations of the invariant distribution of the slow variable S = Xi + 2 X 2 
of system with parameters | |5.2[ |, through computing the null space of the truncated 

generator of the full system (histogram), of the effective generator computed using the CM A 
(solid line) and computed using the QSSA (dashed line). 



QSSA 

CMA 

TTn 

Relative P difference 

6.347 X 10-^ 

1.796 X 10-2 

- 

LH peak position 

20 

20 

20 

LH peak height 

5.378 X 10-3 

1.591 X 10-2 

1.582 X 10-2 

RH peak position 

309 

295 

295 

RH peak height 

6.192 X 10-2 

4.060 X 10-3 

4.006 X 10-3 


Table 4: Differences in the accuracy of the QSSA and CMA approximations of the invariant 
density of S, with respect to the approximation . 


5.2.3. Comparison of approximation of invariant densities 

One approach to quantifying the accuracy of these two methods of approx¬ 
imating effective generators of the slow variable, is to compare the invariant 
distributions of the two systems with that of the marginalised density of the 
slow variable in the full system. We consider the approximation ttq of the in¬ 
variant density of the full system, truncated to the region O, = {{xi,X 2 ) G 
{0,1,..., 800} X {0,1,..., 1200}, as shown in Figure(a). We can marginalise 
this density to find an approximation of the invariant density of the slow vari¬ 
able, as is shown by the histogram in Figure]^ 

The CMA approximation of the invariant density of the slow variable is 
indistinguishable by eye from the highly accurate approximation computed in 
this manner, as shown in Figure The QSSA approximation, on the other 
hand, incorrectly approximates both the placement and balance of probability 
mass of the two peaks in the distribution. The difference in the quality of these 
approximations is stark. This example is an extreme one, as the parameters 
have been chosen to demonstrate how far apart these two approximations can 
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be, but since the CM A has no additional costs associated with it, the advantages 
of this approach are signihcant. The relative P errors of these two approaches, 
when compared with the approximate density ttq , are given in Table along 
with the position and heights of the two local maxima in the densities. 

The CM A computed the generator on the domain S' G [0, 2000] in around 
55 seconds, and the eigensolver took less than a tenth of a second to hnd the 
null space to approximate the invariant density. This is negligible in comparison 
with the cost of exhaustive stochastic simulation of the full system. 

5.2.4- Conditioned path sampling using effective generators 

Given an approximation of the effective generator of the slow variables, com¬ 
puted using the CMA or the QSSA, we can now employ the methodology of m, 
as summarised in Section o to sample paths conditioned on their endpoints. 
This time, a full eigenvalue decomposition of the matrix Ai — + I was 

computed, so that matrices V and D could be found with V unitary and D 
diagonal, with AA = V~^DV. Then rows of Ad’' = V~^D'^V can be efficiently 
and accurately computed, as required in Q and ([^. 




Figure 5: (o.) 8 trajectories of the slow variable S = Xi + 2 X 2 sampled conditioned on 

5(0) = 20, 5(10) = 195, 5(1) S = {0,1,..., 500}Vt S [0,5] for the system ( |11[ | with 
parameters ( |5.2[ | , using the CMA approximation of the effective generator, (b) The means 
and standard deviations of 100 paths sampled using the QSSA (blue plots) and CMA (red 
plots). 


Figu re [5] presents results using this approach. An effective generator for the 
system ( |liy was computed for the domain Xi + 2 X 2 = S' G 12 = {0,1,..., 500}, 
using both the QSSA and CMA, and then fed into the conditioned path sampling 
algorithm. Figure (a) shows 8 samples of conditioned paths approximated 
using the CMA. Notice that as the transition time between the two favourable 
regions is relatively short compared with the length of the simulation, the time 
of the transition varies greatly between the different trajectories. This indicates 
that we are producing trajectories with a fair reflection of what happens in a 
transition between these regions. Figure [^(b) shows the means and standard 
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deviations of 100 paths sampled for both methods of computing the effective 
generator. The QSSA, which overestimates the value of the second peak in the 
invariant density, has a higher mean than the CMA. This demonstrates again 
that errors in approximating the effective generator has a knock-on affect to 
applications such as conditioned path sampling. 

The effective generator was computed on the domain S' € [0, 500] for the 
path sampling, which took the CMA close to 5 seconds to approximate. The 
calculation of the probabilities in Q, and the full eigenvalue decomposition of 
the generator matrix on this domain, took around 50 seconds. After this, each 
path took around 350 seconds to sample. 

5.3. An Example of the Nested CMA Approach 

In this section, we will illustrate how the nested approach outlined in Section 
l^can be applied. We will consider an example for which we know the invariant 
distribution of the slow variables. This gives us a way of quantifying any errors 
that we incur by applying the nested CMA and QSSA approaches. 


Ri : 

0 ^ Xi, 

i?2 : 

X3 ^ 0, 

i?3 : 

^ X2 

i?4 : 

X2 ^ All 

i?5 : 

X2 ^ W3 

i?6 : 

^3 ^ a:2 


We will consider this system in the following parameter regime: 

kiV = 20, k 2 = l, K = 100, 7 = 10. (17) 

As before, V denotes the volume of the well-mixed thermally-equilibrated reac¬ 
tor. In this regime, there are multiple different time scales on which the reactions 
are occurring. This is demonstrated in Figure where there is a clear gap in 
the frequency of reactions i?i and i ?2 (the slowest), i ?5 and Rq (fast reactions) 
and i ?3 and R 4 (fastest reactions). 

This system was chosen as we are able to, using the results in [52], find the 
exact invariant distribution of the slow variable = Xi + X 2 + X 3 . In this 
instance, it is a Poisson distribution with mean parameter 

fci V 

A = - -( 7 fc 2 -I- 37 K -I- 2k2K) = 64.2. 

k2^K 

5.3.1. QSSA-based analysis 

One method to analyse such a system would be a nested QSSA-based analy¬ 
sis, similar to that which is suggested in m- In this paper the authors consider 
systems with reactions occurring on multiple timescales. If at first we consider 
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Figure 6: Relative occurrences of the reactions Ri-Rq, for the system l |16| l with parameters 
The most frequent reactions are reactions R 3 and R 4 , reactions R^ and Rq are the next 
most frequent, with reactions Ri and R 2 being the least frequent. 


all reactions Rs-Rq to be fast reactions, then by applying the QSSA we are 
interested in finding the invariant distribution of the following fast subsystem: 


Cl 

Xi 

+ A2 

+ A 

Rs 

Xi 

K ^ 

^2, 

R4 

X2 

K ^ 

^1, 

As 

X2 

7 ^ 

^3, 

R& 

A3 

7 ^ 

A2. 


( 18 ) 


Note that the quantity Si = Xi + X 2 + X^ is a conserved quantity with 
respect to these reactions, and as such is the slow variable in this system. This 
is in itself also a system with more than one timescale, and as such, we may 
want to iterate again and apply a second QSSA assumption, based on the fact 
that reactions R 3 and i ?4 are fast reactions in comparison with reactions i ?5 
and Rq. This leads to a second fast subsystem: 


Cl : Xi + X2-\-X^ — Si, 

C2 : Xi + X2 = S2, 

R3 ■■ ^ X2, 

R4 : X2 ^ Xi. 

Note that the quantity S 2 = Ai + A 2 = — A 3 is a conserved quantity with 

respect to these reactions, and as such is the slow variable in this system. At this 
point in m, the authors simulate the system using the Gillespie SSA. We could 
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adopt the approach that we used in Section 5.2 in which we find the invariant 
distribution of the system by constructing its generator and then finding the 
normalised eigenvector corresponding to the null space of that operator. This 
would not be expensive since there are only S 2 different states. However, as 
in Section |5.1[ as this system only contains monomolecular reactions, we can 
exactly find its invariant distribution. In this case, Xi and X 2 follow a binomial 
distribution with mean Xi , X 2 = ^ ■ This can then be used to compute the 
effective rate of reaction R 5 in the first subsystem (18), a 5 {Xi,X 2 ) « 7 X 2 = 
^5. This fast subsystem is then reduced to the following: 


Cl 

C2 

i?5 

i?6 


Xi+X2 + X3 = Si= S2+ X3, 
Xi + X2 = S2, 

~fl2 


S 2 

X 3 


X3, 

S 2 . 


Note that we have completely eliminated the fast variables Xi and X 2 , and 
instead consider the slower variable S 2 = Xi + X 2 , with effective rate for R 3 
given by the analysis above. This system is exactly solvable, and its invariant 
distribution is a gamma distribution with means given by X 3 = ^ and S 2 = 
^1^, found by computing the steady states of the mean field ODEs[22]. This 
in turn can be used to compute the effective rate of reaction i ?2 in the full 
system, where we now lose all of the fast variables Xi , X 2 , X 3 and instead wish 
to understand the dynamics of the slow variable S'! = + X 2 + X 3 , which is 

only altered by reactions i?i and i? 2 - This system is given by the following: 


Ri ■. 0 A 5 i, 

R2 : 0 . 


Here the effective rate for R 2 has been found by using the approximation of the 
effective rate q; 2 («S'i) = ^ 2-^3 = ^S. 

5.3.2. CMA-based analysis 

We will now go through the same procedure, but this time using the con¬ 
strained subsystems instead of the fast subsystems as used in the previous sec¬ 
tion. There are 3 choices for the fast reactions, each involving two out of Xi, 
X 2 and X 3 . Since Xi is the product of a zeroth order reaction, it is preferable 
not to include this as one of the fast variables, and so we pick Fi = [X 2 ,X 3 ]. 
We then construct the constrained subsystem for this choice of slow and fast 


24 




variables: 


Cl : 

Xi 

+ X, 

+ X3 — s 

R2 ■■ 

X3 


Xi, 

Rs ■■ 

Xi 

K, ^ 

X2, 

R 4 : 

X2 

K ^ 

Xi, 

R5 : 

X2 

7 ^ 

X3, 

Re : 

X3 

7 ^ 

X2. 


Note that Ri is removed, since it does not change the fast variables. i ?2 is the 
only other reaction which has changes to its stoichiometry due the constrained 
stoichiometric projector. We have reduced the dimension of the system (due 
to the constraint Xi + X 2 + ^3 = a for some cr G N), but we are still left 
with a multiscale system, which in theory could be computationally intractable 
for us to find the invariant distribution for, through funding the null space of 
its generator. Therefore, we can apply another iteration of the CMA to this 
constrained system. 

Reactions i ?3 and R 4 are the fastest reactions in the system, and therefore 
we pick our next slow variable that we wish to constrain to be S 2 = Xi + A 2 , 
which is invariant with respect to these reactions. Due to the previous constraint 
Si = Xi + X 2 + X 3 , we are only required to define one fast variable for this 
system. Both choices F 2 = Xi,X 2 , are essentially equivalent, and so we pick 
F 2 = Xi. These choices lead us to the following second constrained system: 

Cl 

C2 
i ?2 

i?3 

i?4 



ifX2 > 0, 

otherwise. 


( 20 ) 


Here Vi denotes the stoichiometric vector associated with reaction Ri, i.e. the 
vector which is added to the state X(t) if reaction Ri fires at time t. Notice 
that we now have two separate constraints, and as such reactions R 3 and Rq 
now have zero stoichiometric vectors. Moreover, these constraints lead us to 
a somewhat unphysical reaction for i? 2 - The reactant for this reaction is X 3 , 
but only X 2 and Xi are affected by this altered reaction. In system (191 when 
reaction R 2 fires, we lose one X 3 , and gain Xi. Therefore, both constraints 
within (20) have been violated. In order to reset these constraints, without 
changing the fast variable F = X 3 , we arrive at the stoichiometry presented 
in (20). Note that we add the condition that this reaction can only happen if 
X 2 > 0 , as we cannot have negative numbers of any species. 
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This is a closed system, with a very limited number of different states. There¬ 
fore, it is computationally cheap to construct its generator, and to find that 
generator’s null space. Our aim with this system, is to find the invariant distri¬ 
bution of the fast variable given particular values for the constraints Ci and C 2 - 
This distribution will then allow us to compute the expectation of the reaction 
i ?4 within the constrained system ([^, which is the only reaction which is depen¬ 
dent on the results of the second constrained system (since X 3 = Si — 82 )- Once 
the invariant distribution has been found, this can be used to find the effective 
propensity of reaction R 5 given values of = X 1 +X 2 + X 3 and S 2 = Xi+X 2 . 
In turn, the constrained system (191 can then be solved to find the invariant 


distribution on X 3 given a value of Si. Finally, this leads us to the construction 
of an effective generator for the slow variable S*!. 

Since this final constrained subsystem is aphysical, we cannot use the results 
of |22j to find the invariant distribution, and as such we must approximate them 
through finding the null space of the generator, as we did in Section [5.2| 


5.3.3. Comparison of approximation of invariant densities 
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Figure 7: Approximations of the invariant distribution of the slow variable S = Xi -\-X 2 +Xs 
of system with parameters 0- through marginalisation of the invariant distribution of 
the full system (histogram), of the effective generator computed using the CMA (solid line) 
and computed using the QSSA (dashed line). 



Figure shows the invariant distributions of the slow variables S = Xi + 
X 2 + X 3 computed by marginalising the invariant distribution of the full system, 
and from the CMA and QSSA as outlined above. The distribution computed 
using the CMA is indistinguishable by eye from the true distribution, and has 
a relative error of 5.936 x 10“^^, which can be largely attributed to rounding 
errors and error tolerances in the eigenproblem solvers. The QSSA approxima¬ 
tion, on the other hand, has a significant relative error of 3.739 x 10“^. This 
demonstrates again the substantial improvements in accuracy that we gain in 
using the constrained approach rather than one based on the QSSA. This is 
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delivered at no substantial additional computational effort. As in the previous 
two examples, the highly accurate effective generator approximated using the 
CMA can be used in a host of applications where the full generator could not, 
such as conditioned path sampling. 

The CMA is more expensive in this example than the previous ones, as there 
are a very large number of small eigenvalue problems to solve. This is due to the 
fact that there are reactions of three species occurring on three different time 
scales. The generation of the CMA approximation of the effective generator 
took around 1240 seconds, and the subsequent approximation of the invariant 
distribution of the slow variables took just over half a second. This still pales into 
comparison with the cost of exhaustive stochastic simulation of the system. The 
savings would be even more pronounced in systems with multimodal invariant 
distributions where switches between the modes are rare. 

6. Conclusions 

In this paper, we presented a significant improvement and extension to the 
original constrained multiscale algorithm (CMA). Through constructing and 
finding the null space of the generator of the constrained process, we can find 
its invariant distribution without the need for expensive stochastic simulations. 
The CMA in this format can also be used not just to approximate the param¬ 
eters of an approximate diffusion, but to approximate the rates in an effective 
generator for the slow variables. 

In this paper we have not discussed how the slow and fast variables in these 
systems can be identified. In the simple examples presented, this is relatively 
straightforward. However in general, this is far from the case. If the high 
probability regions in the statespace are known a priori, or possibly identified 
through short simulations of the full system, then it is possible to identify which 
are the fast reactions in the system, and therefore what good candidates for 
the slow variable(s) could be. Other more sophisticated approaches exist, for 
example methods for automated analysis to identify the slow manifold [lU 
126) . One relatively ad hoc approach might be to briefly simulate the full system 
using the Gillespie SSA, which can give a good indication as to which the fast 
reactions are. Good candidates for slow variables are often linear combinations 
of the species who are invariant to the stoichiometry of the fast reactions, as we 
have seen in this paper. If the regions which the system is highly likely to spend 
the majority of its time are known, then looking at the relative values of the 
propensity functions, as we did in Figure]^ (b), can lead to an understanding of 
which reactions are fast and which are slow. 

Through iterative nesting, the CMA can be applied to much more complex 
systems, as it can be applied repeatedly if the resulting constrained system is 
itself multiscale. This makes it a viable approach for a bigger family of (possibly 
biologically relevant) systems. This nested approach breaks up the original task 
of solving an eigenvalue problem for one large matrix per row of the effective 
generator, down into many eigenvalue solves for significantly smaller generators 
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for smaller dimensional problems, making the overall problem computationally 
tractable. 

In the hrst example, we demonstrated that the CMA produces an approx¬ 
imation of the dynamics of the marginalised slow process in the system which 
is exact, at least by the measures that we have applied thus far, in the case of 
systems of monomolecular reactions. Since such systems are well understood, 
we were also able to compare this with the accuracy of the equivalent QSSA- 
based method, which incurred significant errors. We then applied the method 
of Fearnhead and Sherlock[T6] to the approximate effective generators of the 
two approaches, in order to approximately sample conditioned paths of the slow 
variables. This task would be computationally intractable to attempt with the 
full generator for this system. This also demonstrated how the accuracies of the 
two approximations can impact the accuracy of any application for which they 
may be used. 

In the second example, a more complex bistable system was also analysed 
using the CMA, and the invariant distribution of the computed effective gen¬ 
erator was shown to be very close to the best approximation that we could 
make of the invariant distribution of the slow variables, using the null space 
of the original generator with as little truncation as we could sensibly manage 
with our computational resources. This was in stark contrast with the poor ap¬ 
proximation which was computed using the equivalent QSSA-based approach. 
This highlighted again the improvement, at no or little extra cost, of using the 
constrained approach as opposed to the QSSA. 

In the final example, we demonstrated how the constrained approach might 
be applied to a more complex example with multiple timescales. The algo¬ 
rithm can be applied iteratively in order to reduce the constrained subsystems 
themselves into a collection of easily solved one-dimensional problems. When 
comparing the invariant distributions of the approximate processes computed 
using the two approaches, the QSSA once again was incorrectly approximating 
the distribution of the fast variables conditioned on the slow variables, and so 
incurred significant errors. In contrast, the CMA produced an approximation 
to the invariant measure which was accurate up to 12 digits. 

We showed how these effective generators can be used in the sampling of 
paths conditioned on their endpoints. Such an approach could be employed as 
a method to sample missing data within a Gibb’s sampler when attempting to 
find the structure of a network that was observed[TH]. This approach could also 
be used simply to simulate trajectories of the slow variables, in the same vein as 
[Sjor[n]. In this instance, it would only be necessary to compute the column of 
the effective generator corresponding to the current value of the slow variables. 

The constrained approach consistently significantly outperforms approxima¬ 
tions computed using the more standard QSSA-based approach, and at negligi¬ 
ble additional cost. Furthermore, in the limit of large separation of timescales, 
the constrained approach asymptotically approaches the QSSA approximation. 

The computational savings that we make in using the CMA depends on the 
application with which we wish to use the effective generators. Similarly, if we 
wish to approximate the invariant distribution of the slow variables, then the 
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CMA will always be less costly than exhaustive stochastic simulation. This is 
because we are able to directly compute the invariant distribution, whereas in 
the simulation setting, to obtain the same statistics we would be required to 
compute a very long simulations. 

If, on the other hand, we simply wish to use the CMA to compute a tra¬ 
jectory of the slow variables, then the savings will vary, based on the size of 
the chosen domain, and the relative differences in propensity of the fast and 
slow reactions in the relevant regions. If our aim is only to produce one rel¬ 
atively short trajectory, then it is possible that stochastic simulation will be 
more efficient than using the CMA. However this is such a trivial task, that any 
modeller wishing to do so what not consider invoking any approximations such 
as the QSSA or CMA. 

There are many avenues for future work in this direction, not least its appli¬ 
cation to more complex biologically relevant systems. In particular, the treat¬ 
ment of systems where the effective behaviour of the slow variable(s) cannot be 
well approximated by a one-dimensional Markov process need to considered, for 
example systems which exhibit oscillations. Automated detection of appropriate 
fast and slow variables, and statistical tests for the validity of the approxima¬ 
tion for different systems would be hugely beneficial. In the case of constrained 
systems which are deficiency zero and weakly reversible, using the results of [T] 
we can find the invariant distributions without even constructing the generator, 
and this could be a good direction to investigate. 
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